Stationary phase-kink states and dynamical phase transitions controlled by 
surface impedance in THz wave emission from intrinsic Josephson junctions 
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As possible states to characterize THz wave emission from intrinsic Josephson junctions without 
external fields, the McCumber-like state and 7r-phase-kink state have been proposed. In the present 
article it is numerically shown that both states are stationary according to the bias current J and 
surface impedance Z. The McCumber-like state is stable for low J and small Z. For higher J, 
the 7r-phase-kink state accompanied with symmetry breaking along the c axis is stable even for 
Z = 1, though strong emission in the vicinity of cavity resonance points only takes place for larger 
Z. Different emission behaviors for Z = 1 and 10 are precisely compared. The dynamical phase 
diagram for 1 < Z < 10 and the optimal value of Z for the strongest emission are also evaluated. 

PACS numbers: 74.50.+r, 85.25.Cp, 74.25.Nf 
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Introduction. Emission from intrinsic Josephson junc- 
tions (IJJs) such as Bi2Sr2CaCu208 (BSCCO) has been 
intensively studied as a candidate of stable source of con- 
tinuous terahertz electromagnetic wave. Although emis- 
sion controlled by flow of Josephson vortices^ had been 
investigated numericallyr^^ experimental realization of 
such emission^ had been quite difficult. Recently more 
evident emission from IJJs without external fields was 
reported experimentally^ and two types of states were 
proposed theoretically in order to explain the emission. 
One is the McCumber-like state (Ohmic and transla- 
tional invariant along the c axis)^ and another is the 
7r-phase-kink state (non-Ohmic with nontrivial symme- 
try breaking along the c axis) 4 These proposals share 
the basic framework jisl and difference in results is due 
to choice of the boundary condition. In the present arti- 
cle we start from a simplified version^ of the dynamical 
boundary condition^ used in Ref. 7, and gradually in- 
troduce effect of the surface impedance.— 

Model and formulation. Differently from emission in 
external fields, dimensional reduction along the field can- 
not be justified anymore in zero-field emission. Actually, 
the above two research groups have already generalized 
their results to three dimensions ! 11 : 12 Nevertheless, such 
generalization did not resolve the discrepancy between 
them. Then, in order to clarify the origin of this discrep- 
ancy, we concentrate on the two-dimensional modeling 
assuming uniform solutions along the y axis £2. That is, 
we solve the following differential equations, 3 
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where the subscript "I + 1,Z" denotes quantities in the 
insulating layer between the Z-th and (I + l)-th super- 
conducting layers, and the operator A*- 2 ) is defined in 
A( 2 )AVu = X l+2tl+1 - 2X l+u + The electric 

field and gauge-invariant phase difference controlled by 
the dc bias current J are basic quantities, and the mag- 
netic field is obtained from d x iif>i+i,i = (1 — C,A^)B[ +l { . 



In these formulas the following scaled quantities are used, 

x' = x/X c , t' = u p t, J' = J/J c , (3) 
E't+1,1 = ivcKPJc)) E? +1>1 , B' = (2nX c d/4>o) B, (4) 
C = Ay M, a = e^ 2 /M, (3 = ^7 c a c X c /(e c c), (5) 
= c/ (v / 4 A c) , Jc = o /(2^oA 2 d), (6) 

with the penetration depths A {, = 0.4/im and A c = 
200/^m, thickness of superconducting and insulating lay- 
ers s = 3A and d = 12A, respectively, Debye length /i = 
0.6A, dielectric constant of the junction e' c — e c /e = 10 
with permittivity of the junction e c , plasma frequency 
u) p , conductivity a c , flux quantum </>o, and critical cur- 
rent J c , following the material parameters of BSCCO in 
Ref. 3. They give a — 0.1 and (3 = 0.02 is taken here. 

Width of the junction L x — 86/im is comparable to 
those of samples used in experiments^ Since direct sim- 
ulation of several hundreds of layers is difficult, the peri- 
odic boundary condition along the c axis corresponding 
to effectively infinite layers is introduced instead. In ad- 
dition to calculations for the number of layers TV = 4, 
some systems with N = 8 or 12 are analyzed to confirm 
numerical consistency. In place of considering outside of 
IJJs, the dynamical boundary condition on edges^ is in- 
troduced. For infinite and uniform layers, this boundary 
condition is simplified^ as the relation between the dy- 
namical part of scaled boundary fields B[ +1 1 and E[ +1 ; : 
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with the dielectric constant of dielectrics e^. The factor 
z(> 1) is considered^ to appear when the wavelength A 
of emitted electromagnetic wave is much longer than the 
thickness of junctions L z , which usually holds in exper- 
iments, as z ss X/L z . Effect of impedance mismatch on 
the edges is also included in Z. When surface fields are 
not uniform along the c axis, Z depends on wave number 
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FIG. 1: Voltage dependence of emission intensity for Z = 1. 
Regions (a)~(f) are divided by jumps of intensity (dashed 
lines), and the data points with star symbols correspond to 
the voltages at which Figs. 2 and 3 are drawn. 



and frequency in the exact expression! 10 i 13 i 14 Width of 
the sample L x is divided into 80 numerical grids. 

In numerical evaluation of stationary states, procedure 
for parameter sweep is essential. One possible approach 
is to start from a zero-current state and gradually vary 
currents^ similarly to experiments. However, it may not 
be suitable for study on stationary states. Simulated time 
scale is much shorter than that in experiments, and acci- 
dental trap by metastable condition may be held during 
gradual change of current. Then, we start from random 
configurations for each current and check the consistency 
of results. Continuous data obtained from independent 
initial conditions strongly suggest that the results are 
stationary ones. Since each data point is independent, 
convergence of calculations can be easily checked by addi- 
tional simulations. In order to obtain convergent results 
from random initial configurations precise algorithm with 
automatically optimized time steps is essential, and the 
RADAU5 ODE solver— is utilized for this purpose. 

Even if this precise algorithm is used, self-consistent 
evaluation of (E[ +1 ; ) is still difficult. We first fix the 
static value of the surface electric field in order to ob- 
tain a stationary state under this constraint, and then 
determine (E' l+1 ; ) self-consistently from this stationary 
state. When the determined average value and amplitude 
of the surface electric field are denoted as E m and AE, 
respectively, strong emission state can be obtained from 
an initial static value E- m i satisfying E- m i < E av — AE. 
This condition looks consistent with the fact that strong 
emission state can only be observed in current-increasing 
process during gradual variance of current £ 

Numerical results for Z — 1. This condition is typi- 
cally realized for z — 1 and e' c — e' d , namely no impedance 
mismatch in an infinite system. Voltage dependence 
of intensity, namely the strength of Poynting vector, is 
shown in Fig.[T] Intensity takes maximum at the onset of 
emission, and there exist several abrupt jumps which cor- 
respond to change of modes of electromagnetic waves in 
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FIG. 2: A series of snapshots of the dynamical part of electric 
fields in a layer of IJJ for some typical values of voltages (per 
layer) for Z = 1: (bl) 0.37mV, (b2) 0.90mV, (c) 1.29mV, (d) 
1.99mV, (e) 2.70mV, and (f) 3.66mV. 
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FIG. 3: Gauge-invariant phase differences in all insurating 
layers corresponding to the voltages for Figs. 2(c), 2(d) and 
2(e). 7r/2- and 37r/2-phase kinks are stabilized in Fig. 3(e). 



IJJs. The range of voltage shown in Fig.Q]is divided into 
6 regions (a) to (f), and snapshots of electric field in each 
region are displayed in Figs.^bl) to[2]^f), respectively. 

The retrapping region (a) is characterized by vanish- 
ing voltage for finite bias current. In the region (b), spa- 
tial dependence of electric fields is small and in-phase 
motion occurs in all layers, which resembles the McCum- 
ber stated This McCumber-like state was automatically 
chosen when translational invariance along the c axis is 
assumed^ A series of snapshots of electric fields for the 
lowest and highest voltages in this region are shown in 
Figs.EJbl) and[2]Jh2), respectively. They show that spa- 
tial dependence of electric fields increases as the voltage 
increases. In the region (c), the state with a 7r-phasc kink 
is favored. This state accompanies symmetry breaking of 
phases along the c axis (see Fig. OUc)). Such a state was 
observed for a large and complex value of or for spa- 
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FIG. 4: Voltage dependence of emission intensity for Z = 
10. Intensity is diverging toward the cavity resonance points 
V — l.lAn [mV/layer] (dashed lines) with the integer n to 
specify cavity modes. Other than the retrapping region (a) 
are divided by these voltages, and the data points with star 
symbols correspond to the voltages at which Fig. 5 is drawn. 
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FIG. 5: A series of snapshots of the dynamical part of electric 
fields in a layer of IJJ at intensity peaks for Z = 10: (b) 
2.28mV (n = 2) and (c) 3.43mV (n = 3) with the integer n 
to specify cavity modes. 



tially inhomogeneous previously. In the regions (d) 
and (f), similar symmetry-breaking states with increas- 
ing number of 7r-phase kinks are stable (see Fig. HJd): 
the 2-kink case with period 4). The region (e) is rather 
special, where two non-7r-phase kinks are stabilized. For 
N = 4, three layers have +7r/2- and — 7r/2-phasc kinks 
and one layer has +37r/2- and — 37r/2-ones (see Fig. 3(e)). 
Structure of non-7r-phase kinks depends on N. For ex- 
ample, for N — 6 four layers with ±27r/3-phasc kinks 
and two layers with =F47r/3-ones. Note that it is not the 
purpose of the present study to check TV-dependence of 
emission states in this region in detail. 

Numerical results for Z = 10. Voltage dependence of 
intensity shown in Fig. 0] is diverging in the vicinity of 
the cavity resonance points (dashed lines), which is quite 
different from the behavior for Z = 1 (see Fig. [I}. The 
McCumber-like state never appears, and standing-wave- 
like behavior (see Fig. [5]) occurs near the cavity resonance 
points. Origin of such behavior can be understood from 
the J-V curve shown in Fig. [6] For Z — 1 (shown in the 
inset), the J-V curve (solid line) almost coincides with 
the Ohm's law (dashed line), which means that most of 
input current changes to Joule heat. While for Z = 10, 



FIG. 6: J-V curve for Z = 1 (inset) and Z = 10. Dashed 
lines represent the Ohm's law for the normal current. 



this curve apparently goes away from the Ohm's law in 
the vicinity of the cavity resonance points, where the 
voltage almost saturates when the current increases and 
the excess energy is emitted as electromagnetic waves. 

Numerical results for other Z . Emission behaviors for 
Z = 1 and 10 are quite different, and it is interesting what 
happens for intermediate values of Z. Then, the dynam- 
ical phase diagram in the Z-J plane is given in Fig. [7] 
Apparently, the McCumber-like state (denoted by M) is 
stable only for a limited parameter region: Z < 5 and 
0.05 < J/J c < 0.12. Boundary of the 2 incommensurate- 
phase- kink state (denoted by I2) is evaluated for N = 4, 
and quantitative change may occur for larger N. Impor- 
tant thing is that this phase is stable only for Z < 2, 
and that most regions of the dynamical phase diagram 
are covered by the 7T-phase-kink states (denoted by K n , 
n: number of kinks). Intensity of emission at the Ki- 
K2 boundary becomes comparable to that at the R-M 
boundary (R: retrapping) for Z«3, and boundaries be- 
tween different 7r-phase-kink states becomes almost inde- 
pendent of Z for Z > 3, where strong emission governed 
by the ac Josephson relation is observed. 

Finally, surface-impedance dependence of emission in- 
tensity for larger Z is investigated. That is, stationary 
emission is optimized for various values of Z, which are 
varied from 3 to 10000. Maximum intensity in each cav- 
ity mode is observed, and all the maximum values are 
plotted versus Z in Fig. [5J The strongest emission is ob- 
served at Z w 50 for n = 1, Z « 80 for n — 2, and 
Z ps 100 for n = 3, respectively. 

Discussions. In Ref. 8, large surface impedance \Z\ = 
1000 was introduced as a consequence of large z due 
to small thickness of the IJJ^ though this surface ef- 
fect may be cancelled^ by penetration of magnetic fields 
from transverse directions, which is neglected in two- 
dimensional modeling. Then, stability of standing-wave- 
like behavior for Z > 3 is important, because emission 
from an infinite IJJ to vacuum (e£. = 10, e' d = 1 and 
z = 1) results in Z — %/TO and satisfies Z > 3. 

Z and n dependence of emission intensity can be ex- 
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FIG. 7: Dynamical phase diagram in the Z-J plane. R, M, 
K„, and I m denote the retrapping, McCumber-like, n-7r-phase 
kink, and m-incommensurate-phase-kink states, respectively. 
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FIG. 8: Surface-impedance dependence of maximum intensity 
in a semi-log scale for first 3 cavity modes: n = 1 (triangles), 
2 (squares), and 3 (circles). 

plained as follows: Electromagnetic wave in IJJs becomes 
closer to standing wave for larger Z , and amplitude of the 
electric field increases. On the other hand, the amplitude 



cannot exceed its static value and saturates as Z further 
increases, and intensity decreases because the magnetic 
field is inversely proportional to Z. As larger value of 
n is taken, static value of the electric field increases and 
the optimal point of emission shifts toward larger Z. 

Although phase difference radically varies from layer to 
layer in the phase-kink states as shown in Fig. [3] surface 
fields are almost uniform for any values of Z. This fact 
shows that the assumption with constant Z in Eq. © is 
a good approximation within the present framework. 

Summary. In the present article terahertz wave emis- 
sion from intrinsic Josephson junctions without external 
fields is investigated numerically, and discrepancy be- 
tween two previous studies based on the two-dimensional 
modeling^ is resolved. The surface impedance Z is 
found out to be essential for characterization of emission 
behavior. For Z = 1, the McCumber-like state is stable 
for low dc bias currents as reported in Ref. 7, and this 
emission state is specific to small Z. Even for Z = 1, 7T- 
phase-kink states become stationary for higher currents. 
Since the 7r-phase-kink states require symmetry break- 
ing of phases between insulating layers, it was not ob- 
served in Ref. 7, where in-phase motion was assumed a 
prior. Although large Z is not necessary for stable ir- 
phase-kink states, strong emission from standing-wave- 
like states does not occur for Z = 1, and such emission is 
observed at least for Z > 3. This condition can be satis- 
fied only by the impedance mismatch between IJJs and 
electrodes, namely the emission from an infinite IJJ to 
vacuum gives Z = y/10. Emission behaviors for Z = 10 
arc qualitatively similar to those in Ref. 8 characterized 
by sharp intensity peaks at the cavity resonance points. 
The strongest emission is observed for Z w 50 ~ 100 
slightly depending on the cavity mode n, and such be- 
havior can be explained by saturation of the electric field. 
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